/**************************************************************
Dimitrije Ruzic
March 2018

This program creates a shell for bootstrapping standard errors
***************************************************************/
clear programs
clear all
set more off
set matsize 11000

local data "/FILE PATH GOES HERE/"
local boots "/FILE PATH GOES HERE/"
local output "/FILE PATH GOES HERE/"

/*Path to amoeba*/
sysdir set PERSONAL "/FILE PATH GOES HERE/"
/*Set the number of bootstraps (nrep) here*/
global nrep 1000

/*******************************************************************************
SECTION 1: Simple estimation without bootstrapped SEs
*******************************************************************************/
local S1 = 0

/*******************************************************************************
SECTION 2: (Much longer) estimation with bootstrapped SEs
*******************************************************************************/
local S2 = 0


/*********************************************************
SECTION 1: Simple estimation without bootstrapped SEs
**********************************************************/
if `S1' == 1 {
/*Create dataset to hold summary of estimates*/
!gunzip /FILE PATH GOES HERE/estimation_dataset.dta.gz
use "$boots/estimation_pairs.dta", clear
	local vnames "bl bk by sigma aL aK"
	foreach name of local vnames {
		gen `name' = .
		}
save "$boots/estimation_pairs_LP_ky.dta", replace
quietly sum estimation_pair
	local pmax = r(max)
	
forvalues num = 1/`pmax' {
	/*Open the estimation dataset*/
	use "$boots/estimation_dataset.dta", clear

	/*Keep only the data needed for this round of estimation*/
	keep if estimation_pair==`num'
	
	matrix A = J(1,3,99)
	
		quietly tsset lbd year
		quietly sum L_elasticity_BLS_5y
		global aLprog = r(mean)
		do "$boots/0_gmm_ky.do"
		
		matrix A[1,1] = 1
		matrix A[1,2] = vcoef[1,1]
		matrix A[1,3] = vcoef[1,2]
		
	/*Find mean and relevant percentiles of estimated coefficients*/
	matrix coln A = iter_num bk by
	svmat A, name(col)
	
	keep iter_num bk by
	keep if iter_num!=.
	keep bk by
	gen bl = $aLprog
	gen estimation_pair = `num'
	
	merge 1:1 estimation_pair using "$boots/estimation_pairs_LP_ky.dta"
		drop _m
		order estimation_pair fk period bl bk by
		sort estimation_pair
		replace sigma = 1/by
		replace aL = bl*sigma/(sigma-1)
		replace aK = bk*sigma/(sigma-1)
	
	save "$boots/estimation_pairs_LP_ky.dta", replace
	
	disp("DONE WITH ITERATION `num'")
	}
}

/*******************************************************************************
SECTION 2: (Much longer) estimation with bootstrapped SEs
*******************************************************************************/
if `S2'==2 {
/*Create dataset to hold summary of bootstrapping*/
use "`boots'estimation_pairs.dta", clear
	forvalues num=1(1)$nrep {
		quietly gen bk_`num' = .
		quietly gen by_`num' = .
		}
save "`boots'estimation_pairs_LP_ky_boots$nrep.dta", replace
quietly sum estimation_pair
	local pmax = r(max)
	
	
forvalues num = 1/`pmax' {
	disp("******************************************")
	disp("Pairing `num'")
	disp("******************************************")
	/*Open the estimation dataset*/
	use "`boots'estimation_dataset.dta", clear

	/*Keep only the data needed for this round of estimation*/
	keep if estimation_pair==`num'
	save "`boots'estimation_b1.dta", replace
	
	/*Draw random samples with replacement (treat each set of establishment 
		observations as an independent block)*/
	local i = 1
	while `i'< $nrep+1 {
		use "`boots'estimation_b1.dta", clear
		quietly {
			/*Create dataset with gid and total observation count*/
			quietly keep if gid[_n] != gid[_n-1]
			quietly sort gid
			quietly keep gid count tt
				tempfile obscount
				save "`obscount'", replace
				
			/*For each of the tt units in the sample, randomly 
				sample an integer between 1 and tt*/
			set seed `i'
			quietly gen kk = int(tt*uniform()+1)
			quietly keep kk
			quietly gen test = 99
			quietly rename kk gid
			quietly sort gid
			
			/*Merge in the observation counts for each of the 
				samples units*/
			merge n:1 gid using "`obscount'"
			drop _m tt
			keep if test!=.
			drop test
			
			/*Expand the dataset to include the correct number of 
				observations per unit*/
			egen fill = fill(1 2/3)
			expand count
			sort fill gid
			by fill: gen smalln = _n
			sort gid smalln
			
			/*Merge in the data from the estimation temp file*/
			quietly merge n:1 gid smalln using "`boots'estimation_b1.dta"
			quietly drop _m
			quietly drop if fill==.
			quietly ren lbd oldlbd
			quietly ren fill lbd
			}
			
		/*Estimate parameters on bootstrapped data and save estimates 
			in matrix A*/
		*disp("******************************************")
		*disp("Pairing `num', Bootstrap replication `i'")
		*disp("******************************************")
		quietly tsset lbd year
		quietly sum L_elasticity_BLS_5y
		global aLprog = r(mean)
		quietly do "`boots'0_gmm_ky.do"
		
		quietly gen bk_`i' = vcoef[1,1]
		quietly gen by_`i' = vcoef[1,2]
		
		quietly keep estimation_pair bk_`i' by_`i'
		quietly keep if _n==1
		
		quietly merge 1:1 estimation_pair using "`boots'estimation_pairs_LP_ky_boots$nrep.dta", nogen
		order estimation_pair fk period
	
		save "`boots'estimation_pairs_LP_ky_boots$nrep.dta", replace
		
		local i = `i'+1
			}
		
		
	}
	
erase "`boots'estimation_b1.dta"
}



